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ABSTRACT 



The main properties of the observed extrasolar plan- 
ets are reviewed with respect to their relevance to 
the formation scenario of planetary systems. Results 
of numerical computations of embedded planets in 
viscously evolving disks are presented. Emphasis is 
given to the accretion and migration process. New 
calculations on inviscid disks are shown. 

The second part of the talk concentrates on reso- 
nant planetary systems. Among the observed extra- 
solar systems there are 3 confirmed cases, Gl 876, 
HD 82943 and 55 Cnc, where the planets orbit their 
central star in a low order mean motion resonance. 
Results of numerical simulations modeling the for- 
mation and evolution of such systems are presented. 

Key words: planet formation, numerical hydrody- 
namics, resonant orbits. 
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Figure 1. Eccentricity vs. distance for the extrasolar 
planets. For comparison the solar system planets and 
the two OGLE objects are also included. 



Since the first discovery of extrasolar planets 
around main sequence stars in 1995, their 
number has risen to over 100 planets to- 
day. For an always up-to-date list see e.g. 
http : //www. obspm.f r/encycl/encycl .html, 
maintained by J. Schneider. Apart from two 
discoveries by t he transit method of OGLE-objects 
(|Konacki et al.ll2003 IpT^izler et al.ll200l . ah others 
have been found by radial velocity measurements. 
Among all systems, there are at least 12 planetary 
systems with two or more planets; a summary 
of their propertie s has recently been given by 
iMarcv et al.l (|2003l) . Longer observations to study 
trends in velocity curves may still increase this 
number. Of all stars investigated, planets have been 
found so far around about 7 %. 

The orbital properties of these planetary systems 
display features quite different from those of our 
own solar sys tem. Some statistical p roperties are 
summarized in|Mazeh <fe Zuckeij (|2nn3)J^antos et al.| 
(|2003t ) and additional data may be obtained from 



http: //www. obspm. f r/encycl/encycl .html. As 

an example we display in Fig^the eccentricity ver- 
sus mass distribution of all plants detected so far. 
While the planets close to the star (hot Jupiters) 
tend to have circular orbits due to tidal circular- 
ization the planets further out have a rather broad 
eccentricity distribution which is essentially indis- 
tinguishable from that of binary stars. Addition- 
ally, the cumulative eccentricity-period distribution 
of stellar binaries and planetary systems are also 
nearly identical for peri ods larger than about 10 days 
(|Mazeh &: Zuckeill2002l) . This may imply a common 
origin of the eccentricities, caused for example by 
dynamical gravitational interaction. 

A summary of the main properties of the exo- 
planets is outlined in the following. The reasons why 
those characteristics present problems to the stan- 
dard planet formation scenari os, is given just below 
the individual points, see also iLissauerl (|2002l) . 



Large Masses (0.2-17 Mj^p) 
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Tidal interaction between the planet and 
the protoplanetary disk results for planets 
larger than about 1 Jupiter mass {Mjup) in 
the cre ation of a density £ap in the disk 
()Lin Paoaloi zou 1980, 1993). In turn, this 
gap will reduce the accretion onto the planet 
and limit the ma ss to about 1 M .j„.n. just as in 
our solar system ()Lissauerlll993l) . 

• Smah Distance (0.04 - 3.3 AU) 

The high disk temperature in the vicinity of 
the star makes it difficult to allow for con- 
densation of solid material. Thus, the cores 
of the protoplanets are expected to have been 
formed beyond the so called snow line at a dis- 
tance of about 3 AU a nd larger from the star 
iSasselov T.eca,Hl200?l) . 

• Large Eccentricities (0.0 - 0.91) 

It is usually assumed that planets form in a pro- 
toplanetary disk, which is essentially in a Kep- 
lerian rotation around the star. Hence, directly 
after formation, planets are expected to orbit on 
circular orbits. 

• Resonant Orbits (2:1, 3:1) 

The condensation proce ss in the disk and 
the subsequent Qlig:ar chic (|Kokubo Ida "19981: 
iThommes et al.ll2003l) growth of planetesimals 
towards larger planets can occur theoretically at 
arbitrary orbits. For stability reasons, only the 
planets cannot be on orbits separated by less 
than about 10 Hill radii. 

• Free Floaters (in Jupiter-mass range) 
Recently, objects in the planetary mass range 
(a few Mj.ir,) but not orbiting stars have bee n 
found (iLucas Rochell20M iB^iar et a1.ll20?ll . 

Although not directly part of a planetary sys- 
tem today, the origin of these objects may be 
related to planet formation, followed by ejec- 
tion processes. It is believed that bodies of a 
few Mjup are to small to have formed by direct 
gravitational collapse. 

All the raised features may have their origin in an 
evolution process of the young planets still embedded 
with the protoplanetary disk from which they formed 
originally. 



2. MODELING EMBEDDED PLANETS 



We shall focus here in the joint evolution of planets 
and disks, and concentrate on the fully non-linear 
hydro dynamical calculations. In recent years sev- 
eral simulations have been performed ([Ki ev 19|^ 
iBrvden et aTlll9 99: 'Lub ow et al.l[T999l: iNdson etldJ 
r2000«: >D^Angelo et al 2003) . In the following we de- 
scribe how these calculation typically proceed. 

All models consider essentially only the late stages 
of the planet formation scenario. It is assumed that 
protoplanets have already been formed for example 
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Figure 2. A typical computational grid-system (^128 x 
128^ used in planet- disk computations 

by core accretion. The subsequent evolution of the 
embedded planets is studied considering the mutual 
gravitational interaction between disk and planets. 
In this contribution we shall focus on the purely hy- 
drodynamic evolution, as elsewhere in this volume 
Richard Nelson describes magnetic effects as well. 

The ingredients of a model for embedded planets con- 
sists of a star (of typical one solar mass, 1 M©), a 
protoplanetary disk {Mdisk = O.OIM©), and a given 
number of planets. 

The disk is modeled by a hydrodynamical evo- 
lution in the field of the star and the plan- 
ets. It is is assumed to be geometrically thin 
with a small vertical thickness H/r « 1. 
Mostly it has been modeled by a two- dimensional 
approximation in r — cp-coordi nates (iKlevI 119991: 
Brvden et an[T99a iLu bow et a1. 111995 [Nelson et al.| 
l2000t iD^Angelo et al.ll200a) . However recently, also 
fully three dimensiona l models (r, cp, i^) have been 
performed (iKIev et a1.l 120(111 iD^Angelo et "all l2003 
IB ate et all l2003l) . Tb capture all relevant effects 
global models of the disk extending typically from 
about 1 to 30 AU or more have to be invoked. Ini- 
tial density profiles typically have power laws for the 
surface density S ex with s between 0.5 and 1.5. 
Nearly all published models have a fixed radial tem- 
perature distribution. Mostly, the aspect ratio H/r 
is held constant in the disk. From there the tempera- 
ture profile follows T(r) ex r~^. In three dimensional 
models the temperature is vertically constant and 
the density p follows a Gaussian. 

For the anomal ous viscosity a Reynolds stress ten- 
sor formulation (|Klevlll995 is used where the kine- 
matic viscosity u is either constant of given by an 
a-prescription u = aCgH^ where a is constant and 
Cs is the local sound speed. From observations, 
values lying between 0.001 and 0.01 are inferred 
for the a-parameter. Full MHD-calculations have 
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Figure 3. Surface density profile for an initially ax- 
isymmetric planet- disk model after 200 hundred or- 
bital orbits of the planet. 



shown that the viscous stress-tensor ansatz may give 
(for sufficiently long time averages) a reasonable ap- 
proximation to the mean flow in a turbulent disk 
(fapaloizoufcNelsoA 2003,). 




Figure 4- Surface density and velocity arrows for a 
one Jupiter mass model. 




The embedded planets are assumed to be point 
masses (using a smoothed potential), and together 
with the star they are treated as classical N-body 
system. The disk also influences the orbits through 
the gravitational torques. This is the desired effect 
to be studied which will cause the orbital evolution 
of the planets. The planets may also accrete mass 
from the surrounding disk. 



The initial conditions of the disk are axisymmetric 
with some given profile for E(r) and T(r) as outlined 
above. Then a Jupiter type planet is placed into this 
disk at a distance of several AU (eg. 5.2), and the 
joint evolution of the (hydrodynamic) disk and the 
planetary (N-body) system is then followed through 
numerical integration. In this review we shall not 
give any details on the numerical issue s, and refer 
the reader to the approp riate literature (|Klevlll99M 
119991: iBrvden et al .11 19991 ) . Just to give some impres- 
sion of how a typical (low-resolution) grid system 
looks like we display in Fig|2l the r — (f grid hav- 
ing a resolution of 128 x 128 gridcells. The radial 
coordinate extends from 1.3 AU up to 20 AU. For 
illustration the star is indicated by central dot and 
the Roche-lobe of a 1 Mjup planet by the solid line 
around the planetary position Xp = — 5.2, = 0. As 
can be seen, at this resolution the Roche-lobe of the 
planet is merely resolved by a few grid-cells. 



Figure 5. Structure of a nested grid centered on 
the planet. Three lev els of grids are shown (from 
iD'Anqelo et aj] IMM))- 

3. VISCOUS LAMINAR DISKS 

The type of modeling outlined in the previous sec- 
tion yields in general smooth density and velocity 
profiles, and we refer to those models as viscous lam- 
inar disk models. In contrast to models which do not 
assume an a priori given viscosity and rather model 
the turbulent flow directly. 



3.1. The global view 

A typical result of such a viscous computation ob- 
tained with a 128 X 280 grid is displayed in Fig 13 
Here, the planet with mass Mp = IMjup and semi- 
major axis ap = 5.2AU is not allowed to move and 
remains on a fixed circular orbit. The disk has 
H/r = 0.01 and z/ = 10~^ in units of a^Qp^ where Qp 
is the orbital Keplerian frequency of the planet. This 
viscosity for the given temperature is equivalent to 
a = 4 • 10~^ at the location of the planet. Clearly 
seen are the major effects an embedded planet has on 
the protoplanetary accretion disk. The gravitational 
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Figure 6. Density Structure of a IMjup on each level of the nested grid system, consisting of 6 grid levels in 
total. 



force of the planet leads to spiral wave pattern in the 
disk. In the present calculation (Fig|2|) there are two 
spirals in the outer disk in the inner disk. The tight- 
ness of the spiral arms depends on the temperature 
(i.e. H/r) of the disk. The smaller the temperature 
the tighter the spirals. The second prominent feature 
is the density gap at the location of the planet. It 
is caused by the deposit of positive (at larger radii) 
and negative (at smaller radii) angular momentum 
in the disk. The spiral waves are corotating in the 
frame of the planet, and hence their pattern speed 
is faster (outside) and lower (inside) than the disk 
material. Dissipation by shocks or viscosity leads to 
the deposit of angular momentum, and pushes mate- 
rial away from the planet. The equilibrium width of 
the gap is determined by the balance of gap-closing 
viscous and pressure forces and gap-opening gravita- 
tional torques. For typical parameter of a protoplan- 
etary disk, a Saturn mass planet will begin to open 
a visible gap. 

The density structure and flow field in the vicinity of 
the planet for such a model with the same physical 
parameter but a resolution of 128 x 440 grid-cells, 
is displayed in Fig^ The solid line represents the 
size of the Roche- lobe of the planet. The velocity ar- 
rows are calculated with respect to the planet in the 
corotating frame. At this resolution the Roche-lobe 
is still barely resolved and the rotation of the mate- 
rial around the planet is visible. It is noticeable that 



even for a one Jupiter mass planet with a very deep 
gap, mass accretion is still possible. Only beyond 
about 5 Mjup the gap gets to wide and additional 
accretion is strongly inhibited. 



3.2. The local view 

To obtain more insight into the fiow near the planet 
and to calculate accurately the torques of the disk 
acting on the planet, a much higher spatial resolu- 
tion is required. As this is necessary only in the 
immediate surrounding of the planet, a number of 
nested-grid and also v ariable grid-size simu l ations 
have been perfo rmed (|D'Angelo et al.1 12002L 120031: 
iBate et ani2003D . As an example we display in FiglSl 
the structure of such a nested grid system with three 
levels. The equations are solved on all levels and 
appropriate boundary conditions are formulated to 
ensure mass, energy and angular momentum conser- 
vation across the grids. This type of grid-system is 
not adaptive, as it is defined in the beginning and 
does not change with time. The planet is placed in 
the center of the finest grid. 

The result for a 2D computation using 6 grids 
is d isplaye d in Fig; El for more details see also 
iD'A ngelo et al.l (|2002t ). The top left base grid has 
a resolution of 128 x 440 and each sub-grid has a size 
of 64 X 64 with a refinement factor of two from level 
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to level. It is noticeable that the spiral arms inside 
the Roche-lobe of the planet are detached from the 
global outer spirals. The two- armed spiral around 
the planet extends deep inside the Roche-lobe and al- 
low for the accretion of material onto the planet. The 
nested-grid calculations have recently been extended 
to three dimensions (3D) and a whole range of plan- 
etary masses have been investigated, starting from 1 
Earth mass to a few Jupiter masses (jP'Angelo et al.l 
I2OO3I) . In the 3D case the strength of the spiral arms 
are weaker and accretion occurs primarily from re- 
gions above and below the midplane of the disk. 



3.3. Accretion and Migration 

These high-resolution numerical computations allow 
for a detailed computation of the torque exerted by 
the disk material onto the planet, and its mass ac- 
cretion rates. Fig gives the inverse 1/tm of the 
migration for three dimensional nested grid calcula- 
tions as a function of the planet mass, given in units 
of the mass of the central star, q = Mp/M^. The 
upper dark solid line represe nts an analyt ical two- 
dimensional linear es timate bv | ^ardl (Il997fl and the 
lower straight line bv lTanaka et al.l (|2002t ). which as- 
sumes a three-dimensional flow and takes the corota- 
tion torques into account. The symbols refer to dif- 
ferent approximations of the potential of the planet. 
It can be seen that for low masses q ^ 10~^ and inter- 
mediate masses q ^ S - 10~^ the numerical results fit 
well to the linear theory. In the intermediate range 
of about q ^ 3 ■ 10~^ the migration rates are about 
an order of magnitude longer ( D'Angelo et al. 200^. 
This effect may be caused by the onset of gap forma- 
tion in the mass range of about 10-15 earth masses. 
Here nonlinear effects begin to set in and modify the 
physics. Thes e resu lts should be compared to those 
obtained bv iiBate et al. ( 2 003 ). 

The consequences of accretion and migration have 
been studied by numerical computations which do 
not hold the planet fixed at some radius but 
rather follow the orbital evolution of the planet 
In Fig[~ 



I we display the radial 
evolution of an embedded protoplanet as it evolves 
with the disk. The typical migrations timescales are 
found to be of the order 10^ yrs, while the accre- 
tion timescale may be slightly smaller. When planets 
reach the vicinity of the central star they can reach 
up to about 4 Mjup. 



4. INVISCID DISKS 




Figure 7. The inverse of the migration rate for differ- 
ent planet masses. The symbols denote different ap- 
proximations (smoothening) for the potential of the 
planet. 



allowed to move and for stability purposes its mass 
is gradually switched on during the first 5 orbits. 
Results are displayed at 4 different times (5, 10, 
15, 32). After only 5 orbits the spiral waves are 
already clearly visible, as they form on a dynamical 
time scale. The gap is just beginning to clear. At 
the edges of the gap high density blobs {vortices) 
are forming which are reduced in number through 
merging. Eventually only one big blob remains. 
Also inside the gap, some detailed structure is 
visible. In comparison, fully viscous simulation 
show these vortices mostly as transient features 
in the beginning of the simulations. They vanish 
later on as a result of the viscosity. If accretion 
and also migration of the planet is considered, these 
moving vortices may create some more complex 
time dependence. However, the details of such 
inviscid models will have to be studied in future 
work. Vortices in accretion disks have sometimes 
been considered to enhance planet formation either 
through triggering a direct gravit ational instabil- 



ity or by trapping particles in it iGodon Livid 
Jf200a); 'de la Fuente Mar cos Barerel (|200lD : 
[Klahr Bodenheimer ( 2003,). 



To investigate the influence of viscosity we studied 
recently also inviscid models with no physical 
viscosity added. Only a numerical bulk vis cosity 
has to be added to ensure numerical stabilitv iKlevI 
However, this has no influence on the 
physical effects of the simulations. In FiglHl we 
display the density structure of such an inviscid 
model for a 1 Mj^p planet. The planet is not 



5. RESONANT SYSTEMS 

In this part we would like to concentrate on plane- 
tary systems which are in resonance. The above men- 
tioned hydrodynamic simulations with single planets 
have been extended to mod els which contain multiple 
planets. It has been shown (|Klevll200nl:lBrvden et al.l 
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Figure 8. The radial position of an initially IMjup 
planet evolving with the disk. In one of the cases the 
mass was held fixed while in the other mass accretion 
from the disk was taken into account. 




Figure 9. Results of an inviscid disk computation 
with an embedded 1 Mjup planet. Results are dis- 
played at four different times (5, 10, 15, and 32) in 
units of the orbital period of the planet. Shown are 
gray-scale plots of the surface density in an (r — cp)- 
coordinate system. 



2^ l^nell^rove et iNelson fc Papain 

2002|) that during the early evolution, when the plan- 
ets are still embedded in the disk, different migration 
speeds may lead to an approach of neighboring plan- 
ets and eventually to resonant capture. More specif- 
ically, the evolution of planetary systems into a 2:1 
resonant configuration was seen in the calculations 
prior to the discovery of any such sys- 
tems. 

In addition to hydrodynamic disk-planet simu- 
lations, many authors have analyzed the evolu- 
tion of multiple-planet systems with N-body meth- 
ods. Each of the known re sonant s y stems have 
been consid e red in detail. iJi et al.l (|2002t ) and 
iLee Peald ()2002aD have modeled the evolution of 
2:1 resonant system GJ 876, while th e 3:1 sy stem 
55 Cnc has been analyzed bv iJi et all (|2003H ) and 
lLee_ fc Peale" (2QQ2b y and the 2 :1 system H P 82943 
by Kjoz dziewski & Macieiewskil (}2QQ1) and iJi et al.l 
()2003al) . Based on orbit integrations, these papers 
confirm that the planets in these systems are in res- 
onance with each other. The dynamics and stability 
of resonant planet ary systems in gene ral has been 
recently studied by iBeauge et al.l (|2002t ) . 

Here we present numerical calculations treating the 
evolution of two planets still embedded in a proto- 
planetary disk. We use both hydrodynamical sim- 
ulations and simplified N-body integrations to fol- 
low the evolution of the system. In the first ap- 
proach, the disk is evolved by solving the full time- 
dependent Navier- Stokes equations simultaneously 
with the evolution of the planets. Here, the motion of 
the planets is determined by the gravitational action 
of both planets, the star, and the disk. In the latter 
approach, we take a simplified approximation and 
perform 3-body (star plus two planets) calculations 
augmented by additional (damping) forces which ap- 
proximately ac count for the gravit ational influence 
of the disk (e.g.lLee fc Peali[2002^. Using both ap- 
proaches, allows a direct comparison of the alterna- 
tive methods, and does enable us to determine the 
damping parameters required for the simpler (and 
much faster) second type of approach. 



5.1. Observations 



The basic orbital parameters of the three known 
systems in mean motion resonance are pre- 
sented in Table ^ The orbital parameters for 
Gl 876 are taken f r om t he dynamical fit of 
iT-^ntrhlin Chamh^ arid for HP 82943 

from iGozdziewski &: MacieiewskJ ((200 ll) . Due to 
the uncertainty in the inclinations of the systems, 
Msini, rather than the exact mass of each planet, is 
listed. By including the mutu al perturbations of the 
planets into their fit of Gl 876,|Laug;hlin <^ ^hamberj 
((20011 ). however, are able to constraint that system's 
inclination to ~ 30° — 50°. 



Two of the systems, Gl 876 and HD 82943, are in 
a nearly exact 2:1 resonance. We note that in both 
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Table 1. The orbital parameters of the three systems 
known to contain a mean motion resonance. P de- 
notes the orbital period, M sin i the mass of the plan- 
ets, a the semi-major axis, e the eccentricity, and w 
the angle of periastron. It should be noted that the 
orbital elements for shorter period planets undergo 
secular time variations. Thus in principle one should 
always state the epoch correspondinq to these osculat- 
ing elements (see e.g. Laughlin & Chambers 2001). 



Name 


P 


Msinz 


a 


e 






[d] 


[Mjup] 


[AU] 




[deg] 


Gl 876 


(2:1) 










c 


30.1 


0.56 


0.13 


0.24 


159 


b 


61.02 


1.89 


0.21 


0.04 


163 


HD 82943 


(2:1) 










b 


221.6 


0.88 


0.73 


0.54 


138 


c 


444.6 


1.63 


1.16 


0.41 


96 


55 Cnc 


(3:1) 










b 


14.65 


0.84 


0.11 


0.02 


99 


c 


44.26 


0.21 


0.24 


0.34 


61 


d 


5360 


4.05 


5.9 


0.16 


201 



cases the outer planet is more massive, in one case 
by a factor of about two (HD 82943) and in the other 
by more than three (Gl 876). The eccentricity of the 
inner (less massive) planet is larger than that of the 
outer one in both systems. For the system Gl 876 the 
alignment of the orbits is such that the two periastrae 
are pointing in nearly the same direction. For the 
system HD 82943 these data have not been clearly 
identified, due to the much longer orbital periods, 
but they do not seem to be very different from each 
other. The third system, 55 Cnc, is actually a triple 
system. Here the inner two planets orbit the star 
very closely and are in a 3:1 resonance, while the 
third, most massive planet orbits at a distance of 
several AU. 



5.2. Modeling resonant planets 

The set of coupled hydrodynamical-N-body models 
presented in this contribution are calculated in the 
sa me manner as the models described above and 
in iKlevI (|l998l [l99 9h for single planets and in iKlevI 
()200n[ ) for multiple planets. The reader is referred to 
those papers for details on the computational aspects 
of the simulations. Other similar models, following 
explicitly the motion of single and multiple planet s 
in disks, have been p resen ted bvlNelson et a 1 (120001) , 
iBryden et al.| and l^nelle^rove ejaU (|2flill|) . 



The initial hydrodynamic structure of the disk, 
which extends radially from rmin to rmax^ is axisym- 
metric with respect to the location of the star, and 
the surface density scales as Hfr) = S pr"-^^^, with 
superimposed initial gaps (|Klevl f2(KK]l) . The initial 
velocity is pure Keplerian rotation {vr = 0, 'U^^ = 
GM^jr^l^^. Here, we assume a fixed temperature 
law with T(r) oc which follows from the assumed 
constant vertical height ajr = 0.10. The kinematic 
viscosity v is parameterized by an a-description v — 
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Figure 10. Overview of the density distribution after 
1250 orbital periods of the inner planet. Higher den- 
sity regions are brighter and lower ones are darker. 
The star lies at the center of the white inner re- 
gion bounded by rmin = 1 AU. The location of the 
two planets is indicated by the white dots, and their 
Roche-lobes are also drawn. Clearly seen are the ir- 
regular spiral wakes generated by the planets. Regular 
intertwined spiral arms are seen only outside of the 
second planet. 



aCsH^ with the sound speed Cg = Hv^p/r. 

In case of the simplifi ed N-body co mputations, we 
follow the approach of iLee Pealel |2002a.) and pa- 
rameterize the imposed damping by specifying a 
damping rate a/ a for the outer planet only, as this 
is the only planet still in contact with the disk. This 
value may depend on time, taking the disk dissipa- 
tion into account. The eccentricity damping time 
scale is a fixed fraction K of the semi-major axis 
damping, i.e. e/e = Ka/a. Thus, in this case only 
a 3 body system (star and two planets) is followed 
by numerical integration, while the effects of the disk 
are taking into account simply through imposed ad- 
ditional forces refiecting the damping of a and e. 



5.3. Results 



At the start of the simulations both planets are 
placed into an axisymmetric disk at ai = 4 and 
a2 = 10 AU, where the density is initialized with 
partially opened gaps superimposed on an otherwise 
smooth radial density profile. Upon starting the evo- 
lution the two main effects are: 



a) Because of the accretion of gas onto the two 
planets the radial region in between them is de- 
pleted in mass and finally cleared. This phase 
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typically takes only a few hundred orbital peri- 
ods. At the same time the region interior to the 
inner planet loses material due to accretion onto 
the central star. Thus, after an initial transient 
phase we typically expect the configuration of 
two planets orbiting within an in ner cavity o f 
the disk, as seen in FigCSl see also lKlevI ()200fll) . 

b) After initialization, the planets quickly (within 
a few orbital periods) excite non-axisymmetric 
disturbances, viz. the spiral waves, in the disk. 
In contrast to the single planet case these are 
not stationary in time, because there is no pre- 
ferred rotating frame. The gravitational torques 
exerted on the two planets by those density per- 
turbations induce a migration process for the 
planets. 
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Figure 11. The semi-major axis (a), eccentricity (e) 
and position angle of the orbital periastron (w) for 
the two planets versus time. In this example, the 
planets have fixed masses of 3 and 5 Mjup, and are 
placed initially at 4 cind 10 AU, respectively. The 
inner planet is denoted by the black line, the outer by 
the light gray line. The dotted reference line (labeled 
3:1), indicates the location of the 3:1 resonance with 
respect to the inner planet. 

Now, the planets' relative positions within the cavity 
have a distinct infiuence on their subsequent evolu- 
tion. As a consequence of the clearing process, the 
inner planet is no longer surrounded by any disk ma- 
terial and thus cannot grow any further in mass. In 
addition, it cannot migrate anymore, because there 
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Figure 12. The evolution of the difference in peri- 
astron fAtu = VD2 — '^i) vs time for the two planets. 



is no torque-exciting material left in its vicinity. All 
the material of the outer disk is still available, on 
the other hand, to exert negative (Lindblad) torques 
on the outer planet. Hence, in the initial phase of 
the computations we observe an inwardly migrating 
outer planet and a stalled inner planet with a con- 
stant semi-major axis (see the first 5000 yrs in the 
top panel of Fig [TT|) . 

This decrease in separation between the planets in- 
creases their gravitational interactions. Once the ra- 
tio of the planets' orbital periods has reached a ratio 
of two integers, i.e. they are close to a mean motion 
resonance, resonant capture of the inner planet by 
the outer one may ensue. Whether or not this does 
actually happen depends on the physical conditions 
in the disk (e.g. viscosity) and the orbital param- 
eters of the planets. If the migration speed is too 
large, for example, there may not be enough time 
to excite the resonance, and the outer planet wi ll 
continue migrating inward fe.g. lHaghighipoinll 19991) . 
Also, if the initial eccentricities are too small, then 
there may be no capture, particularly for second- 
order resonances such as the 3:1 resonance (see e.g. 
iMurrav Dermotti 119991) . The typical time evo- 
lution of the semi-major axis (a), eccentricity (e) 
and direction of the periastron (tu) are displayed in 
Figtm The planets were initialized with zero eccen- 
tricities at distances of 4 and 10 AU in a disk with 
partially cleared gaps. 

In the beginning, after the inner gap has completely 
cleared, only the outer planet migrates inward, and 
the eccentricities of both planets remain relatively 
smah, (< 0.02). After about 5000 yrs the outer 
planet has reached a semi-major axis with an orbital 
period three times that of the inner planet. The pe- 
riodic gravitational forcing leads to the capture of 
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the inner planet into a 3:1 resonance with the outer 
one. This is indicated by the dotted reference hne 
(labeled 3:1) in the top panel of Fig tTTl which marks 
the location of the 3:1 resonance with respect to the 
inner planet. 

We summarize the following important features of 
the evolution after resonant capture: 



a) In the course of the subsequent evolution, the 
outer planet, which is still driven inward by the 
outer disk material, forces the inner planet to 
also migrate inwards. Both planets migrate in- 
ward simultaneously, always retaining their res- 
onant configuration. Consequently, the migra- 
tion speed of the outer planet slows down, and 
their radial separation declines. 

b) Upon resonant capture the eccentricities of 
both planets grow initially very fast before set- 
tling into an oscillatory quasi-static state which 
changes slowly on a secular time scale. This 
slow increase of the eccentricities on the longer 
time scale is caused by the growing gravitational 
forces between the planets, due to the decreas- 
ing radial distance of the two planets on their 
inward migration process. 

c) The ellipses/periastrae of the planets rotate at a 
constant, retrograde angular speed zu. Coupled 
together by the resonance, the apsidal preces- 
sion rate zu for both planets is identical, which 
can be inferred from the parallel lines in the bot- 
tom panel of FigtTTl The orientation of the or- 
bits is phase-locked with a constant separation 
Atz7 = W2 — zui. The rotation period of the el- 
lipses (apsidal lines) is slightly longer than the 
oscillation period of the eccentricities. 

The capture into resonance and the subsequent li- 
bration of the orbits is illustrated further in FigtT2l 
As suggested in Fig^J (bottom panel) the perias- 
trae begin to align upon capture in the 3:1 reso- 
nance. Initially, during the phase when the eccen- 
tricities are still rising (between 5 and 20 thousand 
yrs), the difference of the periastrae settles interme- 
diately to Aw ~ 180°. Then, upon saturation after 
about 20,000 yrs, the system re-adjusts and eventu- 
ally establishes itself at Aw ~ 107°, with a libration 
amplitude of about 7°. This behaviour can be un- 
derstood by an analysis of the interaction Hamilto- 
nian for resonant systems (^Beauge et al. 2002). By 
minimizing the interaction energy, the equilibrium 
values for Aw (and 6i, see below) can be obtained 
as a function of the mass and eccentricity of the two 
planets. 



5.4. Damped N-body computations 

To illustrate the strengths of a simplified N-body 
computation we display in Fig^J the eccentricity 
evolution of the full hydro and the damped 3-body 
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Figure 13. The evolution of the eccentricities of 
the inner planet (upper curves) and the outer planet 
(lower curves) for the full hydro model (black curves 
covering the whole time range) and a simplified 
damped N-body model (light gray), using a damping 
constant of K = 2.5. 



case. Despite some differences which we attribute to 
the unknown eccentricity damping mechanism, the 
overall agreement is reasonable. For a given semi- 
major axis damping rate, the final values obtained 
for ei and 62 at larger times depend on the initial 
values for the eccentricities and the amount of ec- 
centricity damping. In this case we used an initial 
e{to) = 0.01 for both planets and an eccentricity 
damping factor K = 2.5, i.e. a slightly shorter damp- 
ing time scale for eccentricity as for semi- major axis. 
For all models we find that the eccentricity damp- 
ing rate is of the same order as the semi-major axis 
dampin g, i.e. K = 0(1). This finding is in con- 
trast to lLee Pe ale (2002a) who determined a much 
shorter eccentricity damping time, based on models 
for GJ 876. 



5.5. Resonant angles 

To investigate the mass dependence of the resonant 
capture process we ran a sequence of damped N-body 
models with identical initial conditions but different 
masses. The results in terms of the resonant angle 
Aw = W2 — w\ are displayed in FigtT^ For small 
masses < 2Mjup we find capture occurs only into 
the 2:1 resonance. The low masses do not allow for 
a sufficiently strong interaction at the 3:1 resonance, 
and the outer planet migrates through that point. 
All 2:1 resonances settle into the complete symmetric 
configuration Aw = 61 = (FigEJ. For the 2:1 
resonance Oi = 2A2 — Ai — tui, where denotes the 
mean longitudes of the planets. 

For the 3:1 resonances we see anti-symmetric config- 
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Figure 14- Results of a sequence of damped N- 
body simulations. Plotted is the difference of the 
periastrae Azu = W2 — voi of the two planets af- 
ter capture into resonance versus planet mass, where 
rrip = m2 = mi. The diamonds indicate capture in 
2:1 resonance, while the plus signs are for 3:1 reso- 
nance. The other parameters are fixed, as described 
in the text. The horizontal lines indicate values ofO, 
±100° and ±180°. 

urations with anti-aligned periastrae, |Atz7| = 180° 
and 6i = for masses around rrip = 2Mjup, (when 
the first 3:1 cases begin to occur), while for ah 
larger masses we find preferentially the previous non- 
symmetric configurations, where Aw ^ 110^ and 
I Oil ^ 145°. Here, for the case of 3:1 resonance, 
©1 = 3A2 — Ai — 2u7i. As mentioned above, this be- 
haviour may be an indication of a bifurcation in the 
st ability propert i es of resonant systems, as claimed 
bv lBeauge et all (|2002D . With additional damped N- 
body calculations we also find that systems entering 
a 5:2 resonance (not shown) exhibit a typical anti- 
symmetric behavior \Aw\ = 180° and 61 =0°. 



6. DISCUSSION 



Through high resolution numerical computations of 
planets still embedded in protoplanetary disks it is 
possible to calculate the accretion migration rate 
of planets for the whole mass range from a few 
earth masses to over 1 Jupiter mass. The inferred 
timescales can be significantly shorter than the typ- 
ical lifetimes of the disk. 

Thus, we arrive at the problem of stopping the mi- 
gration close to the stars, as well as for Jupiter in 
our own solar systems. In the case of the close-in 
planets it is believed that tidal interaction with the 
central star possibly through magnetic fields may be 
responsible for braking the planet. In the case of 
the solar system, early dissipation of the solar neb- 



ula must have occurred to prevent the planets from 
moving further in. 

Even though a lot of progress has been achieved dur- 
ing the last years, several problems will have to be 
addressed in the future. Among those are inclined or- 
bits which may lead to disk warping (see /3 Pic disk), 
the inclusion of radiative transport to study obser- 
vational effects, planets in binaries to investigate the 
constraints on the formation. The study of turbu- 
lent (MHD) and radiative disks has just begun and 
possibly will lead in the future to a change of some 
aspects of the migration and evolution scenario of 
young planets. 
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